-- ---------------------------------------------------------------------
-- Configuration file
-- Concentrations: mol/L
-- Length:  microns
-- Surface: microns^2
-- Volume:  microns^3
-- ---------------------------------------------------------------------

local pow = math.pow
local cos = math.cos
local sin = math.sin
local PI  = math.pi;
local exp = math.exp

-- ---------------------------------------------------------------------
-- global parameters
-- ---------------------------------------------------------------------
dt_max = 0.01; -- force correct chemistry every dt_max at most

-- ---------------------------------------------------------------------
-- go from 0 to 1 
-- ---------------------------------------------------------------------
function smooth(X)
	return 3*X^2 - 2*X^3
end

-- ---------------------------------------------------------------------
-- make a continuous bump
-- starting a t_start
-- rising during time_up
-- lasting during time_keep
-- decreasing during time_down
-- ---------------------------------------------------------------------
function bump(t,t_start,time_up,time_keep,time_down)
	local t_up = t_start + time_up;
	local t_pl = t_up    + time_keep;
	local t_dn = t_pl    + time_down;
	
	if( t < t_start ) then
		return 0;
	end
	
	if( (t >= t_start) and (t<t_up) ) then
		return smooth( (t-t_start)/time_up );
	end 
	
	if( (t>=t_up) and (t<t_pl) ) then
		return 1;
	end
	
	if( (t>=t_pl) and (t<t_dn) ) then
		return 1-smooth( (t-t_pl) / time_down );
	end
	
	return 0;
end

-- ---------------------------------------------------------------------
-- make a continuous bump
-- starting a t_start
-- rising during time_up
-- ---------------------------------------------------------------------
function plateau(t,t_start,time_up)
	local t_up = t_start + time_up;
	
	if( t < t_start ) then
		return 0;
	end
	
	if( (t >= t_start) and (t<t_up) ) then
		return smooth( (t-t_start)/time_up );
	end 
	
	return 1
end




function relax(t,down_rate)
	local X = down_rate * t;
	return exp( -X );
end

-- smooth up, keep plateau, exponential decay
function xbump(t,t_start, time_up, time_keep, down_rate )
	local t_up = t_start + time_up;
	local t_pl = t_up    + time_keep;
	
	if( t < t_start ) then
		return 0;
	end
	
	if( (t >= t_start) and (t<t_up) ) then
		return smooth( (t-t_start)/time_up );
	end 
	
	if( (t>=t_up) and (t<t_pl) ) then
		return 1;
	end
	
	-- t >= t_pl
	return relax(t-t_pl,down_rate);
end

-- ---------------------------------------------------------------------
-- Geometric parameters
-- ---------------------------------------------------------------------
--Radius  = 15; -- microns
--Surface = 4*PI*Radius^2;
--Volume  = (4*PI/3.0)*Radius^3;

Volume,Surface = geom_data(12.5,5,5);


-- -----------------------------------------------------------------------------
-- Electric Parameters 
-- -----------------------------------------------------------------------------
Cm = 10 * 1.0e-15; -- Farad/ micron^2
Em = -60e-3;       -- Resting potential in V

-- ---------------------------------------------------------------------
-- Reduced Permeabilities, argument is zeta=FV/RT
-- ---------------------------------------------------------------------

-- Surface * Permeability in micron^3 / s
function SP_K( x )
	return ((3.70714279) * exp( (0.4310085034) * x ) + (1.507274188));
end

-- Surface * Permeability in micron^3 / s
function SP_Cl(x)
	return (1.194531961) * exp( (0.05658342092) * x );
end

-- Surface * Permeability in micron^3 / s

function SP_Na(x)
	return (1.360058333) * exp( (0.05223441992) * x );
end

function SP_H(zeta)
	return 0;
end


-- ---------------------------------------------------------------------
-- thermodynamical constants
-- ---------------------------------------------------------------------

Kw      = 1.0e-14;
K_Henry = 29.41;   -- atm/(mol/L)

K1      = 4.45e-7;  
K2      = 5.6e-11;


function P_CO2(t)
	local  P0    = 40.0;
	 return P0/760.0;
end

function kappa(t)
        local  CO2_aq = P_CO2(t) / K_Henry;
        return CO2_aq * K1;
end

pKY    = 6.2

-- ---------------------------------------------------------------------
-- species with algebraic charge and permeabilites
-- ---------------------------------------------------------------------
library=
{ 
--  name, charge, surfacic permeability (function name | nil )
	{ "H+",     1, "SP_H"   },
	{ "HO-",   -1,  nil     },
	{ "HCO3-", -1,  nil     },
	{ "CO32-", -2,  nil     },
	{ "HY",     0,  nil     },
	{ "Y-",    -1,  nil     }, 
	{ "Na+",    1,  "SP_Na" },
	{ "K+",     1,  "SP_K"  },
	{ "Cl-",   -1,  "SP_Cl" },
	{ "HA",     0,  nil     },
	{ "A-",    -1,  nil     },
}


-- ---------------------------------------------------------------------
-- chemical system
-- ---------------------------------------------------------------------
chemsys = 
{ 
--  name, constant(t), { nu_1, species_1 }, { nu_2, species_2}... 
-- a single { nil } => stop
	{ "water",  Kw ,               {  1, "H+" }, { 1, "HO-"}   },
	{ "bicarb", "kappa",           {  1, "H+" }, { 1, "HCO3-"} },
	{ "carb",   K2,                {  1, "H+" }, { 1, "CO32-"}, { -1, "HCO3-" } },
	{ "buffer",  10^(-pKY),        { -1, "HY" }, { 1, "Y-" }, {1, "H+" } },
	{ "lactate", 10^(-3.86),       { -1, "HA" }, { 1, "A-" }, {1, "H+" } }
}

-- ---------------------------------------------------------------------
-- initializing conditions inside
-- ---------------------------------------------------------------------
init_in =
{
	{ 10^(-7.2),   { 1, "H+" }  },
	{ 00e-3,       { 1, "HY"  }, { 1, "Y-" } },
	{ 0.0000,      { 1, "HA"  }, { 1, "A-" } },
	{ 140e-3,      { 1, "K+"  } },
	{ 10e-3,       { 1, "Na+" } },
	{ 20e-3,       { 1, "Cl-" } }
}

-- ---------------------------------------------------------------------
-- initializing conditions outside
-- ---------------------------------------------------------------------
init_out =
{
	{ 10^(-7.4),    { 1, "H+"  } },
	{ 0,            { 1, "HY"  }, { 1, "Y-" } },
	{ 0.0000,       { 1, "HA"  }, { 1, "A-" } },
	{ 4e-3,         { 1, "K+"  } },
	{ 140e-3,       { 1, "Na+" } },
	{ 100e-3,       { 1, "Cl-" } }
}

-- ---------------------------------------------------------------------
-- enzymatic constants 
-- ---------------------------------------------------------------------

function default_efficiency(t)
	return 1;
end

-- AE
KAE = 10e-3;
efficiencyAE = default_efficiency;

-- NHE
L0 = 1000;
Kr = 1.8e-8;
Kt = 3.6e-6;
efficiencyNHE = default_efficiency;

-- Lactate parameters
LaVm = 1e-3/60.0;  -- vmax
LaKm = 30e-3;      -- Km


function lactateProduction(t)
	return 0;
end


-- NaK-ATPase
-- -----------------------------------------------------------------------------
-- JNak in A/micron^2, argument is zeta=FV/RT
-- -----------------------------------------------------------------------------

function _Shape( X )
  return 0.5 * ( 1.0 + math.tanh(X) )
end


-- fitted from article in Ampers
function INaK_exp( x )
  return 231.318e-12 * _Shape( 0.388429*( x - (-3.29312) ) ); 
end 


Capa_exp = 207e-12;       -- from article 207 pF
Surf_exp = Capa_exp / Cm; -- in micron^2

print( '-- Surf_exp = ' .. Surf_exp )

-- electrical flux in Ampers per micron^2
function JNaK( zeta_m )
	return INaK_exp(zeta_m)/Surf_exp;
end

efficiencyNaK = default_efficiency;




runtype = os.getenv("RUNTYPE");


if ( runtype == "ischemia" ) then
	Ta      = 15 * 60;
    Tchange = 1 * 60;
	function lactateProduction(t)
		local prod_rate = (1e-3)/60.0; -- 1mM/min
		--if (t>=0) and (t <= Ta) then
		--	return prod_rate;
		--end 
		return prod_rate * bump(t,0,Tchange,Ta-2*Tchange,Tchange);
	end
	
	function lactateBump(t)
		return bump(t,0,Ta/2,0,Ta/2);
	end
	
end


if ( runtype == "acidosis" ) then


function P_CO2(t)
	local P0 = 40/760.0;
	local P1 = 48/760.0;

	local grow_time = 1 * 60;
	local keep_time = 0 * 60;
	local down_time = 5 * 60;
	
	local start_time = 1 * 60;
	return P0 + (P1 - P0 ) * bump(t,start_time,grow_time,keep_time,down_time);
end

end

if ( runtype == "acidosis2" ) then


function P_CO2(t)
	local P0 = 40/760.0;
	local P1 = 48/760.0;

	local grow_time = 1 * 60;
	local keep_time = 0 * 60;
	--local down_time = 5 * 60;
	
	local start_time = 1 * 60;
	return P0 + (P1 - P0 ) * plateau(t,start_time,grow_time);
end

end

